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Abstract. We develop a two dimensional, self-consistent, compressible fluid model 
to study evolution of Alfvenic modes in partially ionized astrophysical and space 
plasmas. The partially ionized plasma consists mainly of electrons, ions and signif- 
icant neutral atoms. The nonlinear interactions amongst these species take place 
predominantly through direct collision or charge exchange processes. Our model 
uniquely describe the interaction processes between two distinctly evolving fluids. 
In our model, the electrons and ions are described by a single fluid compressible 
magnetohydrodynamic (MHD) model and are coupled self-consistently to the neu- 
tral fluid via compressible hydrodynamic equations. Both plasma and neutral fluids 
are treated with different energy equations that adequately enable us to monitor 
non adiabatic and thermal energy exchange processes between these two distinct 
fluids. Based on our self-consistent model, we find that the propagation speed of 
Alfvenic modes in space and astrophysical plasma is slowed down because these 
waves are damped predominantly due to direct collisions with the neutral atoms. 
Consequently, energy transfer takes place between plasma and neutral fluids. We 
describe the mode coupling processes that lead to the energy transfer between the 
plasma and neutral and corresponding spectral features. 



1. Introduction 

Alfven waves are ubiquitous in laboratory (Gigliotti et al 2009), space (Gekelman 
2004, Shaikh et al 2009, Shaikh & Shukla 2009a, Shaikh & Shukla 2009b, Shaikh 
& Shukla 2009c, Shaikh & Shukla 2009d) and astrophysical plasmas (Balsara 1996, 
Shukla et al 2006, Shukla & Stenflo 2000). These waves are essentially electromag- 
netic fluctuations in magnetized plasmas that propagate predominantly along an 
ambient or guide constant magnetic field. In the context of partially ionized space 
and astrophysical plasmas, these waves interact with the neutral gas and govern 
numerous properties. For instance, Kulsrud & Pearce (1969) noted that the interac- 
tion of a neutral gas and plasma can damp Alfven waves. Neutrals interacting with 
plasma via a relative drag process results in ambipolar diffusion. Ambipolar diffu- 
sion plays a crucial role in the dynamical evolution of the near solar atmosphere, 
interstellar medium, and molecular clouds and star formation. Oishi & Mac Low 
(2006) investigated the role of ambipolar diffusion to set a characteristic mass scale 
in molecular clouds and find that substantial structure persists below the ambipolar 
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diffusion scale because of the propagation of compressive slow mode MHD waves 
at smaller scales. Leake et al (2005) showed that the lower chromosphere contains 
neutral atoms, the existence of which greatly increases the efficiency of wave damp- 
ing due to collisional friction momentum transfer. They find that Alfven waves with 
frequencies above 0.6Hz are completely damped and frequencies below 0.01 Hz are 
unaffected. They undertook a quantitative comparative study of the efficiency of 
the role of (ion-neutral) collisional friction, viscous and thermal conductivity mech- 
anisms in damping MHD waves in different parts of the solar atmosphere. It was 
pointed out by the authors that a correct description of MHD wave damping re- 
quires the consideration of all energy dissipation mechanisms through the inclusion 
of the appropriate terms in the generalized Ohms law, the momentum, energy and 
induction equations. Padoan et al (2000) calculated frictional heating by ion-neutral 
(or ambipolar) drift in turbulent magnetized molecular clouds and showed that the 
ambipolar heating rate per unit volume depends on field strength for constant rms 
Mach number of the flow, and on the Alfvenic Mach number. 

The neutral gas not only collisionally interacts with plasma, but it also charge 
exchanges with plasma protons in many space plasma applications (Zank 1999, 
Shaikh & Zank 2008). For example, in the local interstellar medium (LISM), the 
low density plasma and neutral H gas are coupled primarily through the process 
of charge exchange. On sufficiently large temporal and spatial scales, a partially 
ionized plasma is typically regarded as equilibrated; this is the case for the LISM. 

On smaller scales, the edge region of a tokamak (a donut shaped toroidal experi- 
mental device designed to achieve thermonuclear fusion reaction in an extremely hot 
plasma) is partially ionized. In a tokmak, the neutral particles result from effects 
such as gas puffing, impurity injection, recombination, charge exchange, and possi- 
bly neutral beam injection processes. The presence of neutrals can potentially alter 
the dynamics of zonal flows and cross-field diffusivity (Singh et al 2004, Groebner 
et al 2001). 

It is worth mentioning that the charge exchange process in its simpler (and 
leading order) form can be treated like a friction or viscous drag term in the fluid 
momentum equation, describing the relative difference in the ion and neutral fluid 
velocities. The drag imparted in this manner by a collision between ion and neutral 
also causes ambipolar diffusion, a mechanism used to describe the Alfven wave 
damping by cosmic rays (Kulsrud & Pearce 1969) and also discussed by Oishi & 
Mac Low (2006) in the context of molecular clouds. 

In this paper we focusc on understanding the propagation characteristic of Alfven 
waves by considering the collisional and charge exchange interactions. Our model 
includes both the interactions simultaneously. In Section 2, we discuss the equa- 
tions of a coupled plasma-neutral model, their validity, the underlying assumptions 
and the normalizations. Section 3 describes the results of our nonlinear, coupled, 
self-consistent fluid simulations. We find that propagation characteristics of Alfven 
waves depends critically on collision and charge exchange processes. Finally, a sum- 
mary is presented in section 4. 

2. MHD Model Equations 

We assume that fluctuations in the plasma and neutral fluids are isotropic, homo- 
geneous, thermally equilibrated and turbulent. A mean or constant magnetic field 
is present. Local mean flows may subsequently be generated by self-consistcntly 
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excited nonlinear instabilities. The boundary conditions are periodic, essentially a 
box of plasma. Most of these assumptions are appropriate to realistic space and 
astrophysical turbulent flows. They allow us to use MHD and hydrodynamic de- 
scriptions for the plasma and the neutral components respectively. In the context 
of LISM and outer heliospheric plasmas, the plasma and neutral fluid remain close 
to thermal equilibirium and behave as Maxwellian fluids. Our model simulates the 
plasma-neutral fluid that is coupled via collisions and charge exchange in space (and 
astrophysical) plasmas. The fluid model describing nonlinear turbulent processes in 
the interstellar medium, in the presence of charge exchange, can be cast into plasma 
density (p p ), velocity (U p ), magnetic field (B), pressure (P p ) components according 
to the conservative form 



dF 



£ + V • Q p = Qp,„, 



(2.1) 
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The above set of plasma equations is supplimcnted by V • B = and is coupled self- 
consistently to the neutral density (p n ), velocity (V„) and pressure (P n ) through 
a set of hydrodynamic fluid equations, 
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Equations (2.1) to (2.2) form an entirely self-consistent description of the coupled 
plasma-neutral turbulent fluid. 

Several points are worth noting. The charge-exchange momentum sources in the 
plasma and the neutral fluids, i.e. Eqs. (2.1) and (2.2), are described respectively 
by terms Qm(U p , V„, p p , p n , T n , T p ) and Qm(V„, U p , p p , p„, T n , T p ). A swapping 
of the plasma and the neutral fluid velocities in this representation corresponds, 
for instance, to momentum changes (i.e. gain or loss) in the plasma fluid as a re- 
sult of charge exchange with the neutral atoms (i.e. Qm(U p , V„, p p , p n , T n , T p ) in 
Eq. (2.1)). Similarly, momentum change in the neutral fluid by virtue of charge 
exchange with the plasma ions is indicated by Qm(V„, U p , p p , p n , T n , T p ) in Eq. 
(2.2). (U p — V„) = — F rlyP is ion-neutral collision force. This force 

is self-consistently calculated in our model. In the absence of charge exchange in- 



4 



D. Shaikh 



teractions, the plasma and the neutral fluid are de-coupled trivially and behave 
as ideal fluids. While the charge- exchange interactions modify the momentum and 
the energy of plasma and the neutral fluids, they conserve density in both the flu- 
ids (since we neglect photoionization and recombination). Nonetheless, the volume 
integrated energy and the density of the entire coupled system will remain con- 
served in a statistical manner. The conservation processes can however be altered 
dramatically in the presence of any external forces. These can include large-scale 
random driving of turbulence due to any external forces or instabilities, supernova 
explosions, stellar winds, etc. Finally, the magnetic field evolution is governed by 
the usual induction equation, i.e. Eq. (2.1), that obeys the frozen-in-field theorem 
unless some nonlinear dissipativc mechanism introduces small-scale damping. 

The underlying coupled fluid model can be non-dimcnsionalizcd straightforwardly 
using a typical scale-length (4>), density (po) and velocity (vo). The normalized 
plasma density, velocity, energy and the magnetic field are respectively; p p = 
p p /po,\J p = \]p/vo,P p = P p / pqVq,~B — B/wo-y/po. The corresponding neutral 
fluid quantities are p n — p n /po,U„ = U n /vo,P n = P n /poVQ. The momentum 
and the energy charge-exchange terms, in the normalized form, are respectively 
Qm = Qm?-o/po v o , Qe = Qeio/po v o- The non-dimensional temporal and spatial 
length-scales are t — ivo/4), x = x/£o- Note that we have removed bars from the 
set of normalized coupled model equations (2.1) & (2.2). The charge-exchange cross- 
section parameter (a) , which does not appear directly in the above set of equations 
(see the subsequent section for more detail), is normalized as a = no^oc, where the 
factor no A) has dimension of (area) -1 . By defining n ,£ through a ce = l/n £ = 
k^ e , we see that there exists a charge exchange mode (k ce ) associated with the cou- 
pled plasma-neutral turbulent system. For a characteristic density, this corresponds 
physically to an area defined by the charge exchange mode being equal to (mpf) 2 . 
Thus the larger the area, the higher is the probability of charge exchange between 
plasma ions and neutral atoms. Therefore, the probability that charge exchange can 
directly modify those modes satisfying k < k ce is high compared to modes satisfying 
k > k ce . Since the charge exchange length-scales are much smaller than the turbu- 
lent correlation scales, this further allows many turbulent interactions amongst the 
nonlinear turbulent modes before they undergo at least one charge exchanges. An 
exact quantitive form of sources due to charge exchange in our model is taken from 
Shaikh & Zank (2008) [see also appendix A]. 

3. Alfven Wave Dynamics 

We have developed a two-dimensional (2D) nonlinear fluid code to numerically in- 
tegrate Eqs. (2.1) to (2.2). The 2D simulations are not only computationally less 
expensive (compared to a fully 3D calculation), but they offer significantly higher 
resolution (to compute inertial range turbulence spectra) even on moderately-sized 
clusters like our local Beowulf. The spatial discretization in our code uses a dis- 
crete Fourier representation of turbulent fluctuations based on a pseudospectral 
method, while we use a Runge Kutta 4 method for the temporal integration. All 
the fluctuations are initialized isotropically with random phases and amplitudes in 
Fourier space. A mean constant magnetic field B is assumed along the y-direction. 
Our algorithm ensures conservation of total energy and mean fluid density per unit 
time in the absence of charge exchange and external random forcing. Additionally, 
V • B = is satisfied at each time step. Our code is massively parallelized using 
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(a) Alfven waves 



(b) Compressible modes 
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Figure 1. Evolution of Alfvenic and fast/slow modes in a partially ionized plasma. The 
effect of charge exchange and collisions is shown, (a) The propagating Alfven waves damp 
predominantly due to collision (than the charge exchange), (b) Compressive (i.e. fast/slow) 
modes show a less effective tendency towards the damping. The lower curves correspond 
to a higher strength of collisional and charge exchange. The modes and time are measured 
in the normalized units £q and la / vq respectively. 



Message Passing Interface (MPI) libraries to facilitate higher resolution. The ini- 
tial isotropic turbulent spectrum of fluctuations is chosen to be close to k~ 2 with 
random phases in all three directions. The choice of such (or even a flatter than 
-2) spectrum does not influence the dynamical evolution as the final state in our 
simulations progresses towards fully developed turbulence. 

Our goal is to determine the quantative evolution of Alfven waves in a partially 
ionized plasma. For this purpose, we distinguish the Alfvenic and non- Alfvenic, i.e. 
corresponding to the compressional or due to slow and fast magnetosonic modes, 
contributions to the turbulent velocity fluctuations. To identify the distinctive role 
of Alfvenic and fast /slow (or compressional) MHD modes, we introduce diagnostics 
that distinguish the modes corresponding to Alfvenic and slow/fast magnetosonic 
fluctuations. Since the Alfvenic fluctuations are transverse, the propagation wave 
vector is orthogonal to the velocity field fluctuations i.e. k _L U, and the aver- 
age spectral energy contained in these (shear Alfvenic modes) fluctuations can be 
computed as (Shaikh & Zank 2006) 



(M*)> 



'Ekl*xu k p 



E k |u k 



The above relationship leads to a finite spectral contribution from the |kx Uk| char- 
acteristic turbulent Alfvenic modes. On the other hand, fast/slow (i.e. compressive) 
magnetosonic modes propagate longitudinally along the velocity field fluctuations, 
i.e. k || U and thus carry a finite component of energy corresponding only to the 
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Figure 2. The rate of transfer of energy (or energy exchange) between the plasma and 
neutral fluids is shown. Clearly the decay rates of energy in one fluid is compensated by 
the increased rates in the other fluid such that entire energy of the coupled plasma-neutral 
system is conserved provided no external source of dissipation exists in the system. 



ik • Uk part of the velocity field, which can be determined from the following rela- 
tionship 



<*c(t)> 



'Ekl*-u k | 



E k |u k | 2 ■ 

The expression of k c essentially describes the modal energy contained in the non- 
solenoidal component of the MHD turbulent modes. 

The evolution of Alfvenic and compressive modes is shown in Fig (1). We vary 
charge exchange strength and collision parameter in our simulation to examine their 
effects on the propagation of Alfvenic and compressive modes in the partially ionized 
plasmas. Clearly, the two processes operate on different time and length scales and 
are self-consistently modeled in our simulations (see Eqs 1 ans 2). We find that 
charge exchange and collisional interactions jointly damp the Alfven waves. This 
is shown in Fig (la) for 256 2 modes in a two dimensional box of length 2-7T x 27r. 
The other parameters in our simulations are; charge exchange k/k ce ~ 0.1 — 0.01, 
fixed time step dt = 10~ 3 , and collision parameter v ~ 0.1 — 0.001. The background 
constant magnetic field B = 0.5. Our simulations are fully nonlinear because the 
ratio of the mean and fluctuating magnetic fields <5B/Bo ~ 1. It is clear from 
Fig (la) that the increasing strength of charge exchange and collision leads to a 
damping of the Alfven waves. Similarly, the compressive modes are also damped 
by virtue of the increased charge exchange and collision. This is shown in Fig (2a). 
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We further find that collisional interactions are more promising to damp the MHD 
waves in comparison with the charge exchange process. 

It is further interesting to note that rate of change of energy in plasma fluid is 
balanced precisely by rate of increase of energy in the neutral fluid by an equivalent 
amount. This is shown in Fig (2). This is true as long as there are no external 
sources of dissipation in the two fluids. We find from Fig (2) that rate of decay 
of plasma is followed by the rate of increase of energy in the neutral fluid. Thus 
on an average, the net transfer of energy from plasma to the neutral fluid (or vice 
versa) is nearly zero. It then turns out that the collisional and charge exchange 
processes dampen the wave activities, while the energy cascade rates in the two 
fluids continue to remain constant in the absence of external dissipation. The effect 
of a random driving does not change the above results, because the random froces 
isotropically drive the plasma and neutral fluids. 



4. Conclusion 

In this work, we develop a two dimensional self-consistent model of plasma and neu- 
tral fluids that are coupled through charge exchange and collisional interactions. 
We concentrate on understanding the propagation of MHD waves in the partially 
ionized space and/or astrophysical plasmas. One of the most important points to 
emerge from our studies is that charge exchange in combination with collisional in- 
teractions modify the propagation characteristics of Alfven and compressive modes. 
The wave activities are damped. By contrast, the energy transfer rates in the plasma 
and neutral components continue to remain constant. It is important to note that 
the collision and charge exchange processes introduce disparate time and length 
scales in the coupled plasma- neutral fluid system. We find that on scales I ^ £ ce , 
the coupled plasma system evolves differently than the uncoupled system where 
large-scale turbulent fluctuations are strongly correlated with charge-exchange. We 
add finally that our self-consistent model can be useful in studying turbulent dy- 
namics of partially ionized plasma in the magnetosphcrc of Saturn and Jupiter 
where outgassing from moons and Io and Encephalus introduces a neutral gas into 
the plasma (Shaikh & Zank 2008). 
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Appendix A. Charge Exchange Sources 

The charge exchange terms can be obtained from the Boltzmann transport equation 
that describes the evolution of a neutral distribution function /„ = /(x, v x ,v y , v z ,t) 
in a six-dimensional phase space defined respectively by position and velocity vec- 
tors (x, v x , v y , v z ) at each time t. Here we follow Pauls et al. (1995) in computing 
the charge exchange terms from various moments of the Boltzmann equation. The 
Boltzmann equation for the neutral distribution contains a source term propor- 
tional to the proton distribution function f p and a loss term proportional to the 
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neutral distribution function /„, 

df F 

-^+v„-V/„ + --V Vn /„ = P-i (Al) 
at m 

P = f p (x,V x ,Vy,V z ,t) [ f n {x,V x ,Vy,V z ,t)\v n - U p \eT cx (v re l)d 3 V n 

L = f n (x,v x ,v y ,v z ,t) J f p (x,v x ,v y ,v z ,t)\u p - v n \(j cx (v re i)d 3 u p . 

Here f p , u p represent respectively the proton distribution function and veloc- 
ity. CTcx is the charge exchange cross-section (between neutrals and plasma pro- 
tons), m is the mass of particle, and F represents forces acting on the fluid. The 
charge exchange parameter has a logarithmically weak dependence on the rela- 
tive speed (v re i — \u p — v n \) of the neutrals and the protons through ct cx = 
[(2.1 - 0.092 ln(v rei ))l(T 7 cTO] 2 [Fite et al 1962]. This cross-section is valid as long 
as energy docs not exceed leV, which usually is the case in the inner/outer helio- 
sphere. Beyond leV energy, this cross-section yields a higher neutral density. This 
issue is not applicable to our model and hence we will not consider it here. The 
density, momentum, and energy of the thermally equilibrated Maxwellian proton 
and neutral fluids can be computed from Eq. (A 1) by using the zeroth, first and 
second moments J f^d 3 £,, J m^f^d 3 ^ and J m£ 2 /2/{<i 3 £ respectively, where £ = u p 
or v n . Since charge exchange conserves the density of the proton and neutral fluids, 
there are no sources in the corresponding continuity equations. We, therefore, need 
not compute the zeroth moment of the distribution function. Computing directly 
the first moment from Eq. (A 1), we obtain the neutral fluid momentum equation as 
given by Eq. (2.2). The entire rhs of Eq. (A 1) can now be replaced by a momentum 
transfer function Qm (««,%)) which reads 

QmK,h p ) = fi{u p ,v n ) - p,(v n ,u p ), (A 2) 

where Qm and fi, the transfer integral, are vector quantities. The transfer integrals 
describe the charge exchange transfer of momentum from proton to neutral fluid 
and vice versa. The first term on the rhs of Eq. (A 1) can be expressed by 

Ji{u p ,v n ) = f p (x,v x ,v y ,v z ,t)l3(up,v n ) 

where, 

f3(u p ,v n ) = J fn(x,v x ,v v ,v z ,i)\v n - u p \a cx (v re i)d 3 v n . 

Considering a Maxwellian distribution for the neutral atoms, we simplify (3{u p ,v n ) 
as follows, 



p{Up,v n ) = a ex n n V Tn i - H T . 2 ■ 

V 71 

Note that the above expression emerges directly from a straightforward integration 
of sources in the rhs of the Boltzmann Eq. (A 1). To obtain the expression for the 
momentum transferred from proton to neutral (or vice versa), we need to take a 
second moment of the p-{u p , v n ) expression. This is shown in the following, 

fi(u p ,Vn) = mv n Io(u p ,v n ) +m(u p - v n )Ii(u p ,v n ). 
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where Iq and I± are transfer integrals that can be written as follows, 

h(u p ,v n ) = J fp(xL,v x ,v y ,v z ,t)[3(u p ,v n ) d 3 u p ; 

h{u p ,v n ) = J v n fp(x,v x ,Vy,v z ,t)(3(u p ,v n ) o9u p . 

Assuming a Maxwellian distribution for plasma protons and using f3{u p ,v n ) from 
the above expression, we can straightforwardly evaluate the transfer integrals Iq 
and ii (see Pauls et al (1995) for details). We further write the complete form of 
the first term on the rhs of Eq. (4) as follows, 



fj,(u p ,v n ) 



U* 



\ n - (u p - v„) 



In a similar manner, we can evaluate the second term on the rhs of Eq. (A 1), which 
yields the following form, 



n(v n ,Up) = mo- cx n p n n 



^*„,« P U P-( V «-U P ) 



Vii .Up 



where 



U* = U* 



u* 



7T 



4 V " 



4 -V? + AU : 



— Vr 



1/2 



-,1/2 



-V4 



-V, 



AW 
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On substituting these expressions in the momentum transfer function, we obtain 



QjuK.Mp) = ma cx n p n n (v n - u p ) 



U* + 



J- n 



(A3) 



Equation (A3) together with Eq. (2.2) yields the momentum equation for the neu- 
tral gas. Swapping the plasma and neutral fluid velocities yields the corresponding 
source term for the proton fluid momentum equation. The gain or the loss in neutral 
or proton fluid momentum depends upon the charge exchange sources, which de- 
pend upon the relative speeds between neutrals and the protons. The thermal speeds 
of proton and neutral gas in Eq. (A 3) are given respectively by V£ p = 2KBT p /m 
(the factor 2 arises because of thermal equilibration in that it is assumed that 
the temperature of the plasma electrons and protons are nearly identical so that 
T p = T e + T pro ton — 2T p ) and Vf, — KsTn/ra. The corresponding temperatures are 
related to the pressures by P p = 2n p KsTp and P. n — n n KBT n , where n„, T„, n p , T p 
are respectively the neutral and plasma density and the temperature, and Kb is 
the Boltzmann constant. 

The moment, J m£ 2 /2f^d 3 ^, of the Boltzmann Eq. (Al) yields an energy equa- 
tion for the neutral fluid whose rhs contains the charge exchange energy transfer 
function 

Q E (Vn,U p ) = T](Up,V n ) - T](v n ,U p ), 
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where r](u p , v n ), r](v n ,u p ) are the energy transfer (from neutral to proton and vice 
versa) rates. These functions can be computed as follows: 



V 2 



1 3 

T)(u p , v n ) = -mV 2 a ex npn n U* + -mVTa ex npn n AVu p ,v n -ma c ^n p n n w n -(u p -v n ) In 

2 4 °v UptVn 

and 

1 3 

r](v n ,u p ) = -mV 2 an p n n U* + -mVx ' cr ex npn n AV Vn , Up -ma ex n p n n Up-(-v n -u p ) 



V T 2 

_ P 



The total energy transfer from neutral to proton fluid, due to charge exchange, can 
then be written as, 

QE(v n ,u p ) = ^ma ex n p n n U* (V 2 - U%) + jma ex n p n n {VS n AVu p ,v n - V^ p AV VnyUp ) 



-ma cx n p n n 



v„ • (u p - v„) 



V 2 



sv Up , Vr 



- Up • (v„ - Up) 



V T 2 



(A 4) 



and 



*V Up , Vn = l-VS p + ^Vl+AU 2 }^, 

IT ' \)TT 



A similar expression for the energy transfer charge exchange source term of plasma 
energy in Eq. (2.1) can be obtained by exchanging the plasma and neutral fluid 
velocities. In the normalized momentum and energy charge exchange source terms, 
the factor ma cx in Eqs. (A 3) & (A 4) is simply replaced by a. 
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